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Abstract. Random matrix theory is used to assess the significance of weak correlations and is well es- 
tablished for Gaussian statistics. However, many complex systems, with stock markets as a prominent 
example, exhibit statistics with power-law tails, that can be modelled with Levy stable distributions. We 
review comprehensively the derivation of an analytical expression for the spectra of covariance matrices 
approximated by free Levy stable random variables and validate it by Monte Carlo simulation. 
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tribution theory and Monte Carlo studies - 02.70.Uu Applications of Monte Carlo methods 



I 1 Introduction 

, The classical ensembles of random matrices play an im- 
• portant role in the modelling of physical systems, in time 
I series analysis and in other fields. The first notion of a 
, matrix ensemble in statistics was given in the 1920s by 
' Wishart for the purpose of correlation analysis [T] . Physi- 
cists began to be interested in random matrices in the 
1950s, when Wigner presented a model of nuclear energy 
[ levels as eigenvalues of symmetric random matrices W 
. whose elements are random numbers drawn from a Gaus- 
' sian distribution N{0,a'^) [2j, or actually from any sym- 
. metric distribution with a finite second moment |3j, e.g. 
■ equiprobable ±1 random numbers. With increasing ma- 
I trix size the eigenvalue spectrum tends to the semicircle 
law: 

Pw(A) - -^V4a2-A2. (1) 

Wigner's data were based on neutron and proton scatter- 
ing. Other applications of random matrix theory in physics 
include classical and quantum chaos, disordered systems, 
many-body quantum systems, quantum dots, quantum 
chromodynamics, quantum gravity, supersymmetric field 
theory, string theory, etc. In 1998 Guhr et al. wrote a re- 
view on many of these with more than 800 references [4] . 
In 2003 the Journal of Physics A dedicated a special issue 
to random matrix theory [5]. Random matrices are used 
in other fields too, e.g. operations research, for diverse 
problems as bandwith efficiency in wireless communica- 
tion [617] or optimal aircraft boarding [8i9j . In correlation 



analysis the theory of random matrices can be used to 
assess whether weak correlations are significant or just 
noise. The mathematical link between correlation matri- 
ces of time series and random matrices is the Wishart 
matrix ensemble, that, together with the Wigner ensem- 
ble, is one of the standard tools in the theory of random 
matrices. Recent introductions to the latter including nu- 
merical aspects can be found in Refs. [lOlllj . Since the 
1990s econophysicists have employed random matrix the- 
ory for the analysis of correlation in financial time series 
I12I13I14I15I16TT7] . with portfoHo theory [l8lT9] as one 
of the motivations; a particular attention is given to the 
largest eigenvalues of the covariance matrix and the asso- 
ciated eigenvectors, that correspond to the whole market 
and its sectors. Recently, random matrix theory was used 
also for a correlation analysis of macroeconomic time se- 
ries [2^. 

Consider i = 1^...,N stochastic time series Xij ob- 
served at synchronous times tj, j = 0, . . . , T. The data 
can be arranged in a TV x T matrix M of increments 
where each row corresponds to a time 
series and each column to a sampling time. Assuming that 
the average of the increments is zero, the Pearson estima- 
tor for the covariance of two time series i and j is 



Cy = T^^mikTrijk. (2) 

fc=i 
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The covariances of all pairs can be collected in a x 
symmetric matrix 



1_. 
T' 



-MM' 



(3) 



The covariance matrix C is also called Wishart matrix as 
it was studied by him. One is often interested in testing 
the hypothesis that there are no significant correlations. 
This can be done comparing the eigenvalue spectrum of an 
empirical correlation matrix with the spectrum of a refer- 
ence matrix built with synthetic uncorrelated time series. 
If the matrix rows are random walks whose increments are 
independent and identically distributed (iid) normal devi- 
ates with standard deviation a, the spectrum describing 
the above null hypothesis in the limit for N, T oo with 
m = N/T is given analytically by the Marcenko-Pastur 
law m\ : 



PC (A) = 



v/(A+-A)(A-A_) 
27rCT^mA 

2 . 



(4) 



A± = 0-2 [i±y^y 



This result has been rediscovered a few times |11I22I23] . 
Indeed, for a sufficiently large matrix the exact distribu- 
tion of its elements becomes less and less relevant, and the 
Marcenko-Pastur law can be obtained for iid increments 
drawn from any symmetric distribution with a finite sec- 
ond moment . This effect was evident also in Wigner's 
studies of matrices whose elements are binary random 
variables assuming the values ±1 with equal probability. 
In both the Wigner and Wishart ensembles the spectra 
of large matrices converge to that of an infinite matrix 
(respectively the semicircle law and the Marcenko-Pastur 
law) as a consequence of a generalised central limit theo- 
rem. 

A practical use of Eq. ^ is that if the empirical spec- 
trum of data shows significant differences from the the- 
oretical curve, then it may be justified to reject the null 
hypothesis of no true correlations. The details of the lat- 
ter are then a separate issue. In principle it is possibile 
to test not only correlation, but also any kind of suit- 
able assumption leading to a given shape of the expected 
spectrum, both theoretically and numerically. Depending 
on the specific case one chooses a suitable null hypoth- 
esis. For example, if the considered time series are the 
log-prices of traded stocks, in a first approximation it is 
reasonable to test the absence of true correlation with nor- 
mally distributed log-returns |12|13|24] . Another powerful 
approach requiring less knowledge of the distribution of 
the increments is a bootstrap scheme that consists in re- 
sampling the covariance matrix after random permuta- 
tions of the empirical time series. Since the reshuffling of 
the rows of M destroys any possible correlation, an ab- 
sence of correlation among the original time series requires 
that the eigenvalue spectrum of C does not change. 

So far, the result given by Eq. (|4]) lies within clas- 
sical random matrix theory and requires iid matrix ele- 
ments with finite moments. In this work we are concerned 
with the Wishart-Levy ensemble as a natural extension of 



the Wishart-Gaussian ensemble treated by the Marcenko- 
Pastur theory. The situation becomes more complicated 
if the elements of M are distributed with power-law tails, 
as happens in numerous physical, biological and economic 
data [21]. Stock markets as well as many other complex 
systems exhibit a dynamics that results in power-law tailed 
statistics. The Marcenko-Pastur theory is not valid any 
more when the second moment is not finite, and the cor- 
responding spectral densities cannot be obtained from a 
simple extension of Gaussian random matrix theory. As a 
consequence of the central limit theorem for scale- free pro- 
cesses the distribution of many of the above phenomena 
is usually assumed to be a symmetric Levy a-stable dis- 
tribution, whose pdf is given most suitably as the inverse 
Fourier (cosine) transform of its characteristic function: 



La{x) — T f 

1 



-1 

k 



7^ 



(x) (5) 

f'^'^)" cos(xfc) dfc. 



The second and higher moments of La (x) diverge for a < 
2, and for a < 1 even the first moment does not exist. 
If a = 2 Eq. ([S]) gives a Gaussian with standard devia- 
tion a = V2j. However, we shall see that the functional 
representation of this distribution is not required in the 
derivation of the spectrum. 

A matrix whose elements are iid samples from a stable 
density is called a Levy matrix. A symmetric Levy matrix 
is called a Wigner-Levy matrix. A symmetric matrix C 
built from a Levy matrix M according to the equation 



C = 



1 



2^2 /a 



MM"^ 



(6) 



is called a Wishart-Levy matrix. Notice that the normali- 
sation factor has been generalised with respect to Eq. ([3]) 
to take into account Levy a-stable statistics. Sampling the 
elements from the probability density function 



fxix)^N'^^LaiN^^"x), 



(7) 



the limiting spectrum becomes independent of the matrix 
size N [25]. It turns out that the spectra of these matri- 
ces have no longer a finite support as in the semicircle 
and Marcenko-Pastur laws and are dominated by the be- 
haviour of the power-law tail of La{x). 

It was proposed to use the theory of free probabil- 
ity with its convenient machinery leading to analytic re- 
sults that could be obtained otherwise only by means of a 
painful use of combinatorics. A free Levy stable random 
matrix has a spectrum belonging to the class of free sta- 
ble laws. The contemporary physical and mathematical 
literature on correlation matrix analysis with power-law 
tailed uncorrelated noise is very active also in the context 
of free probability. Limiting the list to physics journals, the 
reader can consuh Refs. ^ 27 28 29 30I31I32I33I34I35I36I37I38I39] . 
For a review of free probability theory see Ref. [40]. The 
Marcenko-Pastur spectrum can be obtained as a special 
case of this more general theory. 
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Our aim in this paper is to review comprehensively the 
analytic derivation of the spectral density of free stable 
Wishart-Levy random matrices already solved by Burda 
et al. [26127 29 30 31 32 33|34|35j and, as a further step, 
to validate numerically the analytic result by Monte Carlo 
simulation. The rest of this paper is organised as follows. 
Sec. [2] introduces the mathematical background of free 
probability theory, whose objects are elements of an alge- 
bra, usually an operator algebra, and may enjoy the prop- 
erty of freeness. Sec. |3] explains free stability and presents 
an approximation for the Wishart-Levy covariance matrix 
of time series using free stable random variables. An ex- 
planation of free stability is provided too. Sec. [4] derives 
in detail a transcendental equation, due to Burda et al., 
whose solution gives the spectral density for the approx- 
imated covariance matrix. Sec. [5] shows numerically the 
validity of this equation comparing analytical and Monte 
Carlo results. A summary and an appendix with computer 
code conclude the paper. 



2 Mathematical background 

A symmetric NxN matrix X has real eigenvalues Xi, . . . ,Xn 
The spectral density of X can be written as 



1 ^ 



(8) 



where it is assumed that the weight of each eigenvalue is 
the same and each eigenvalue is counted as many times as 
its multiplicity. The resolvent matrix mj is defined as 



Gx(^) = (^i-X)-i, zeC, 



(9) 



where 1 is the NxN identity matrix. The Green function 
is defined as 

Gx(2)-^trGx(z), (10) 

where the trace tr of a square matrix is defined as the sum 
of its diagonal elements. If X is a random matrix, the 
above definition is generalised including an expectation 
operator E: 

Gx(z) = ^E[trGx(z)]. (11) 

The Green function contains the same information as the 
eigenvalues and the eigenvalue density of X [T3]. The Green 
function can be written in terms of the eigenvalues of X: 



1 ^ 1 



N ^ z-X, 

2—1 



(12) 



This is a special case of the definition through the Cauchy 
transform of a generic spectral density: 



By using the following representation of Dirac's (5-function, 

-^^Pv(1)^i7tS{x), (14) 

where PV denotes the principal value, the spectral density 
can be obtained from the Green function: 

px(A) = lim -Im[Gx(A-ie)]. (15) 

e^0+ TT 

This means that the eigenvalues follow from the disconti- 
nuities of Gx(2) on the real axis. 

Non-commutativity of matrices and, in general, of op- 
erators makes it difficult to extend standard probability 
theory to matrix as well as operators spaces. Among possi- 
ble extensions of probability theory to operator spaces the 
so-called free probability theory has the advantage that 
many results can be deduced from well-known theorems 
on analytic functions [34] . 

In order to explain the framework of free probabil- 
ity, let us start from conventional classical probability. 
A probability space (12, JF, P) is a measure space, where 
]7 is the sample space, is a cr-algebra on f2, and P : 

— > [0, 1] e M is a non-negative measure on sets in 
J- obeying Kolmogorov's axioms; ui G f2 is called an el- 
emetary event, A G is called an event. A random vari- 
able X : J7 — > R is a measurable function that maps el- 
ements from the sample space to the real numbers, and 
thus elements from to a Borel cr-algebra S on R. The 
probability distribution of X with respect to P is described 
by a measure fxx on (M, S) defined as the image measure 
of P: fj,x{B) = P[X~^{B)], where B is any Borel set and 
X~^{B) C !F is the counter-image of B. The cumulative 
distribution function of X is Fx (x) = fix {X < x). The ex- 
pectation value for any bounded Borel function g : R — > R 
is 

E[g{X)]= f g{x)fix{dx)= f g{x)dFx{x). (16) 

If Fx{s) is differentiable, the probability density function 
(pdf) of X is fx{x) = dFx{x)/dx. 

This construction can be extended to non-commutative 
variables, e.g. matrices or more in general operators. Let 
A denote a unital algebra over a field F, i.e. a vector space 
equipped with a bilinear product o : ^ x ^ ^ ^ that has 
an identity element I. A tracial state on is a positive 
linear function t : A ^ ¥ with the properties t(I) — 1 
and t(XY) = r(YX) for every X,Y e ^. The couple 
(A, t) is called a non-commutative probability space. 

For our purposes A — B{H), where B{H) denotes the 
Banach algebra of linear operators on a real separable 
Hilbert space Ti.. This is a *-algebra, as it is equipped with 
an involution (the adjoint operation) X t-^ X* : B{Ti.) — > 
B{TL). Considering a self-adjoint operator X e B{H), it 
is possible to associate a (spectral) distribution to X as 
in classical probability. Thanks to the Riesz representa- 
tion theorem and the Stone- Weierstrass theorem, there is 
a unique measure /.tx on (K, S) satisfying 



f+°° 1 
Gx(^)= / rPx(A)dA. 



(13) 



g{x)fiy^{Ax) = r[5(X)] 



(17) 
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where g : M ^ M is any bounded Borel function [iO^. 
Therefore we say that the distribution of X is described 
by the measure /ix- For our purposes this measure is equal 
to the spectral density px defined in Eq. (|15p . In random 
matrix theory the Wigner semicircle law has the role of the 
Gaussian law in classical probability, and the Marcenko- 
Pastur law corresponds to the law. 

Classically, independence between two random vari- 
ables X and Y can be defined requiring that for any couple 
of bounded Borel functions /, g 

E[(/(X) - E[/(X)])(5(r) - E[giY)])] = 0. (18) 

Analogously, two elements X and Y in a non-commutative 
probability space are defined as free (of freely) indepen- 
dent with respect to r, if for any couple of bounded Borel 
functions /, g 



T[(/(X)-r[/(X)])(g(Y)-r[5(Y)])]=0. 



(19) 



Defining freeness between more than two elements is a 
non-trivial extension 42^ . 

Generally, square N x N random matrices X are non- 
commutative variables with respect to the function t(X) = 
(1/A^) E[tr X] , see Eq. ([11]), but for any given N no pair of 
random matrices is free. Nevertheless two random matri- 
ces X, Y can reach freeness asymptotically if for any inte- 
ger n > and any set of non- negative integers (71 , . . . , 7„) 
and (/3i, ...,/?„) for which in the limit N ^ 00 



r(X''i) = . . . = r(X''") = r(Y'3i) = . . . = r(Y^") = 
we have 



(20) 
(21) 



This means that large random matrices can be good ap- 
proximations of free non-commutative variables. 

Given an operator X S B{Ti), the following functions 
are useful in deriving its spectral distribution /ix: 



1. Moment generating function, defined as 
Mx(z) = zGx(z) - 1. 



(22) 



The name stems from the fact that, if the distribution 
of X has finite moments of order fc, mx,fc = t(X'^), 



Mx(z) = 2^^. 



k=l 



(23) 



This can be seen inserting the sum of the geometric 
series 



00 



(24) 



with q = \/\z\ into Eq. ([T 

Gx(^) = — ^— -px(A)dA (25) 

z{l-\/z) 

+00 -, 00 . 



-E;r/'x(A)dA (26) 

E^Tt/ AVx(A)dA (27) 

00 

E^- (28) 



2. R-transform. In classical probability the pdf of the sum 
of two independent random variables X + Y is equal 
to the convolution of the individual pdfs, i.e. 



fx+vix) = {fx * fY)ix). 



(29) 



The convolution is done conveniently in Fourier space, 
where it becomes a multiplication: the characteristic 
function 



fx+Y{k) - / fx+Y{x)e'''' dx (30) 
Jr 

of X + Y is the product of the characteristic functions 
of X and Y, 



fx+Y{k)^fx{k)fY{kh 



(31) 



and the cumulant generating function of X + y is the 
sum of the cumulant generating functions of X and Y: 

log fx+Y (fc) = log fxik)+ log fY{k). (32) 

The free analogue of the cumulant generating function 
is the i?-transform invented by Voiculescu [40 43 44] as 
part of the functional inverse of the Green function: 



Gx ( i?x(z) + i 



(33) 



The i?-transform for the sum of two free operators is 
the sum of their i?-transforms: 



i?X+Y(2) =-Rx(z) + i?Y(2). 



(34) 



The free analogue of convolution is indicated with the 
symbol ffl: 

A^x+Y = Aix ffl A*Y- (35) 
This is computed through i?x, given the connection 
between the Green function Gx and the spectral distri- 
bution /ix- Other definitions of the i?-transform were 
proposed later. 

Blue function. It is convenient to introduce also an in- 
verse of the Green function Gx(2), called Blue function 
as a pun [45] : 

Gx(Bx(^)) = Sx(Gx(^)) = z. (36) 
The Blue function is related to the i?-transform by 



fc=0 



Rx{z) + -. 
z 



(37) 
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4. S -transform. In the same fashion as the i?-transform 
for the sum, another transform allows to compute the 
spectral distribution of the product of two operators 
from their individual spectral distributions: 



where xx(-z) is defined through 



Xx(2Gx(2:) - 1) = -. 



(38) 



(39) 



For X 7^ Y the ^-transform of the product is the prod- 
uct of the individual S'-transforms: 



(40) 



As the i?-transform allows to compute the free addi- 
tive convolution ffl, the 5'-transform leads to the free 
multiplicative convolution Kl: 

MXY = A*x K A^Y- (41) 

3 Free stable random variables and the 
Wishart-Levy ensemble 

Let P be the matrix projector of size T xT, with iV ones 
in arbitrary positions on the diagonal and all the other 
elements zero, e.g.: 



P = diag(... ,1,1,. ..,0,1,0, 0,1,. ..,1,0,. 



(42) 



Let A be a (large) T x T matrix with a free stable spec- 
tral distribution. This property is the analogue of clas- 
sical stability. The sum of two free non-commutative /i- 
distributed variables results in a new ^-distributed vari- 
able. The Wishart matrix ensemble of size N x N defined 
in Eq. ([3]) can be approximated using the N x T matrix 
M/T^/" obtai ned from PA if only the N non-zero rows 
are considered |26I27I29I30I31I32I33I34135] . Indicating this 
operation with curly braces, the approximation reads 



rp2/a 



MM^ ~ {PA}{A"rp}. 



(43) 



The former equation is justified by very good results, both 
analytic and numeric, in a similar approach for Wigner- 
Levy matrices |35) . 

Once we know the domain of attraction for one spe- 
cific classical stable distribution, we can expect that a 
sum of iid random numbers, e.g. Z = (l/A/'n) 
with some suitable normalisation J\fn, converges to their 
attractor for large n. If Zi are independent elements of 
random matrices, as in Ref . [12] , each of them tends to a 
stable law under matrix addition. However, for free sta- 
bility we must consider random matrices as a whole, and 
a different procedure is needed. A fundamental point is a 
property discussed by Bercovici and Pata 46J, that can 
be summarized as follows. If I?c(/^c) and I?f(/if) are the 
domains of attraction of the stable laws /ic and /if in 



classical and free probability respectively, a distribution 
I' € I?c(Mc) V € Vfi^fif). In other words, if we are able 
to recognise the classical attractor of a distribution i/, 
we also know its free attractor Df . Moreover, one and only 
one free stable distribution corresponds to any set of pa- 
rameter values characterising a classically stable distribu- 
tion. The spectrum of a Wigner-Levy matrix is symmetric 
with the same tail index a of its entries, i.e. it belongs to 
the domain of attraction of a well-recognised classical sta- 
ble law. This means that the sum of sufficiently many free 
non-commutative variables with this spectrum converges 
to a non-commutative variable with a stable distribution. 

Another property discussed in Refs. [40l47l48j can be 
summarised for our purpose as follows. Considering two 
N X N matrices L.; and Lj with i ^ j and two inde- 
pendent random orthogonal N x N matrices and Oj, 
the matrices OiLiOj and OjLjOj are free in the limit 
iV — > oo. These properties together with the observation 
that Li and OiLiOj have the same spectrum justify the 
equation [3S] 



1 ^ 
^^;^^0,L,OT. 



(44) 



This means that a free stable non-commutative variable 
can be approximated adding randomly rotated classical 
Levy random matrices. 

To generate Levy matrices we use the Chambers-Mal- 
lows-Stuck algorithm [49 50 : a random number X drawn 
from the symmetric Levy a-stable pdf, Eq. ([5]), can be 
obtained from two independent uniform random numbers 
U,V G (0, 1) through the transformation 



— log U cos 
cos((l - a)<P) 



sin(Q;^) 
cos<? ' 



(45) 



where = 7r(y- 1/2). For a = 2 Eq. (gS]) reduces to X = 
2jy/— log U sm<P, i.e. the Box-MuUer method for Gaussian 
deviates with standard deviation a = V2^. 

The QR-decomposition of a T x T matrix H with ran- 
dom Gaussian entries yields 



H = OU, 



(46) 



where O is random orthogonal and U is upper (or right) 
triangular. For alternative methods to obtain a random 
orthogonal matrix see Ref. |5lj and references therein. 



4 The analytical spectrum 

The moment generating function of the T xT matrix D = 
APA^ satisfies the transcendental equation |26|27|29|34j 

- exp (^i^ j zM'^'^{z) = (AfD(z) + 1)(Md(z) -t- m), 

(47) 

which can be solved analytically for a few special values of 
1/4, 1/3, 1/2, 2/3, 3/4, 1, 4/3, 3/2, 2; the solution 
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was published for a = 1 [27]. The equation can be solved 
numerically for other values, see the Appendix. Actually, 
we are interested in the spectrum of the approximation of 
C provided by the rhs of Eq. (j43|) , but the Green functions 
of the matrices D and C are related by the equation [31] 



Gd(-z) = ™^ Gci'mz) + 



1 — TO 



(48) 



whence, noticing that m Gci'mz) = Gc{z), 

Mr>{z) = zGoiz)-! = mzGciz)-m = toMc(z). (49) 

In the following we will explain in detail the route that 
leads to Eq. (|T7)l and then to the desired spectral density 

As in classical probability stable laws have an analytic 
form for their Fourier transform, free stable laws have an 
analytic form for their Blue transform |35|42|46|52j : 



B\{z; a) ~ a + bz 



(50) 



The parameter a accounts for a horizontal shift in the 
distribution of the matrix elements and can be set to zero 
without loss of generality. The parameter b depends on the 
distribution; for the symmetric Levy a-stable pdf, Eq. ([5]), 
it has the value [29] 



I = gi7r(a/2-l)^ 



(51) 



As discussed in the previous section, given an index a G 
(0, 2], Ba{z; a) indirectly but precisely defines the attrac- 
tor law for the sum of free variables with a-tailed spectral 
distribution. Since free probability theory is exact only 
in the large size limit T,N oo, N/T = to, the only 
variables that define the model are a and to. 

Rewriting Eq. (|50p with Ga(z) in place of z and using 
Eq. ([36]) yields 



(52) 
(53) 



bGT\z) + G^Hz)^z, 
which is equivalent to 

6Q(z) + zGa(z) + 1 -0, Ga{z)^0. 



In Sec. [2] we established calculation rules with the help 
of which the solution of our specific problem can be put to- 
gether piece by piece. First notice that thanks to Eq. (j40]) . 
if for simplicity from now on we substitute A with its sym- 
metrised counterpart (A -I- A^)/2 so that A = A^, 



'APA 



SaS] 



A^PA 



SaS^ 



AJAP 



s. 



AAP 



s 



A2p. 



(54) 



For the ^-transform of the matrix product A^ we also re- 
quire the Green function. The desired relation is a conse- 
quence of the fact that the spectral measure for free Levy 
a-stable operators in the Wigner ensemble is symmetric 



Pa(A) = Pa(-A) 
Ga(^) = G_a(^). 



(55) 
(56) 



The Green function of A^ can be expressed in terms of 
the Green function of A exploiting the Cauchy transform 
representation and the previous symmetry: 



Ga4z) = 



+ 00 



1 



— oo 



z- A2 
1 



PA (A) dA 
1 



2\/i 



2ViVV^-A V^+A 
(GA(Vi) + G_A(x/i)) 



- ^GA(^/I). 

\/z 



PA(A)dA 



(57) 



The next piece in the composition of the solution is the 
5-transform of the projector P, which requires its Green 
function too. Inserting the spectral density of P, 



pp(A) = mS{X - 1) -f (1 - m)S{X), 



(58) 



into the definition of the Green function of P as a Cauchy 
transform yields 

Gp(z) = J j^pp{X)dX 

[mS{X - 1) + (1 - m)5{X)]dX 



1 



z-X' 
TO 1 — m 



z-1 



z 



(59) 



The moment generating function Mp(z) = zGp{z) — 1 
and the definition of the S'-transform finally give 



Sp{z) 



1 



Z + TO 

Rewriting Eq. (|53p with y/z in place of z, 
bGli^) ~ V^Gi(V^) + 1 = 0, 
and inserting Eq. ([57)1 yields 

6z"/^G^2(z) - zGa2(z) + 1 = 0. 
Observing that from Eq. ([35)1 

1 1 



(60) 

(61) 
(62) 



Xa2(2Ga2(z) - 1) 



, (63) 



Eq. (|62|) becomes 
bx-'^'^Gl 



{-) 


Ga2 




\Xa^J 


XA2 


Kxa^J 



Because from Eq. (|38)) it follows that 

-1 = ^, 

2 / 

Eq. can be simplified to 



Ga2 ( 

XA2 VXa2 



1 = 0. (64) 



(65) 



^Xa?/'G^ ( — ) = z. 



1 



(66) 
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Multiplying both sides by X^"^^/^ yields 



1 

XA2 



Z -a/2 



then subtracting and adding 1, 









( — Ga^ I 


(xaO 




VXa2 







-q/2 



and inserting again Eq. (|65p gives 



which can be written as 



XA2 



2/q 



(67) 



(68) 



(69) 



(70) 



(z + l)2 V6. 

Now, using the definition of the S'-transform and the result 



1 + z 1 

-^A^ = XA2 = 



z(l + z)\h 



11 a 



(71) 



which can be used to write 5d, the 5-transform of the 
Wishart matrix on the rhs of Eq. (|43p is 



'5'pA2 — '5'pS'a2 — 



2/a 



z(m + z) 



(72) 



This result is the starting point for the way back. Re- 
applying the definition of the 5'-transform we can write 



XA2p 



z 1 
^ + 1 (z + l)(z + m) 



and 



Xa2 



ip = (^ + l)(^ + ™)(^) 



-2/q 



(73) 



(74) 



Together with Mx^iz) = zGxiiz) — 1 this allows to sub- 
stitute Xd(-Wd(^)) = 1/^ and AfD(l/XD(^)) = z. Notice 
that we changed the index A^P to D to emphasise our 
goal. So we can finally write 



z= (MD(z)-(-l)(MD(z)+m) 



A/d(z) 



-21 a 



(75) 



Inserting Eq. (j49p yields the corresponding equation for 
C: 



mMc(z)\"^/ 



z = (mMc(z) -f l)(mAfc(2) + m) 
gathering m: 

2 = m2-2/"(Afc(2) + l/m)(Afc(2) + 1) 



(76) 

-2/q 

(77) 



From Eq. (|22p and from the relation between the moment 
generating function and the spectrum we finally obtain 



PC (A) = — Im[Afc(A + ^0-)]. 



(78) 



Inserting h from Eq. ([5T|) and rearranging, Eq. ([75]) takes 
the form anticipated in Eq. (|47p . Returning to the moti- 
vation of the paper, the result described by Eq. ([77]) must 
be considered an approximation of the curve correspond- 
ing to the null hypothesis of absence of correlation in time 
series with fat-tailed increments. 



5 Monte Carlo validation 

It has already been shown numerically that the theory 
works in the Wigner-Levy ensemble [35] . For the Wishart- 
Levy case we produced free Levy stable random matrices 
A of size T X T through Eq. (|44|) : a, N x N principal mi- 
nor of AA^ is a free Wishart-Levy matrix C with the 
desired asymmetry ratio m = N/T < 1. Such a minor 
results from the action of the projectors P in Eq. ((43|) . 
Since a square matrix of size T contains n = \T/N \ non- 
overlapping principal minors of size N < T, this procedure 
can be repeated for the same matrix A with different pro- 
jectors Pi, where i — l,...,n labels the projector that 
selects the rows from {i — 1)N -I- 1 to iN. Especially if m 
is small, it is computationally favourable to follow closely 
Eq. (gSl) by first building an N x T matrix = {P^A} 
made of N rows out of A, and then forming the product 
Cj = MiMj. The eigenvalues of Ci are accumulated in 
a histogram that gives the spectrum. This procedure is 
repeated producing enough matrices Ci until the desired 
statistical accuracy is reached. All plots in Fig.[l]have been 
made using an equal number of eigenvalues for the sake 
of comparability. Free stable laws as defined by the Blue 
function in Eq. (|50p and the empirical spectra have differ- 
ent normalisations. For the purpose of a comparison as in 
Fig.IU this is corrected dividing M by a factor r(l-|-Q;)^/", 
that can be obtained comparing the asymptotic behaviour 
of the two spectra. The Appendix gives the code for the 
calculation of the spectral density by Monte Carlo as just 
described. 

This procedure implements the definition of the Wishart 
covariance matrix based on a real random rectangular 
data matrix M. In this paper free probability theory has 
been used to provide an analytic equation for the spec- 
trum of a Wishart matrix with the simplifying assump- 
tion that the matrix A on the right hand side of Eq. (|43)) 
is symmetric. Therefore M may contain symmetric ele- 
ments too, which is not necessary in the definition of the 
Wishart ensemble. However, it is possible to see that this 
does not affect the properties of MM^. In other words, 
the symmetrisation introduced for simplicity in the an- 
alytic derivation does not change the original numerical 
problem by introducing correlations. Actually, our Monte 
Carlo scheme does not symmetrise the matrix A obtained 
from Eq. (|44p and matches the analytic spectrum. 
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Fig. 1. Spectral densities from the numerical solution of the analytic equation (solid lines) and from Monte Carlo simulation 
(stairs). In each case the dimension of C is A'^ = 200, the number of addends in Eq. (|44|) is i? = 20, and the number of sampled 
eigenvalues is 5* = 36 000. 



6 Summary 

We have explained the justification as well as the math- 
ematical basis with which free probability theory enters 
random matrix theory, in particular in the context of the 
Wishart matrix ensemble. Since the derivation of the an- 
alytic solution for the spectra of free stable random ma- 
trices has not been published in a self-contained way yet 
|26I27I29I34I55] . we recollected it in detail. Then we vali- 
dated numerically with Monte Carlo calculations the ana- 
lytic prediction of the eigenvalue spectrum for free stable 
Wishart-Levy matrices. Overall we find an excellent con- 
sistency between theory and simulation. 
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Appendix: Computer codes 

The numerical solution of Eqs. ((77H781) was computed with 
Mathematica 6.0 in almost one fine: 

a = 3/2; 
m = 1/3; 
width = 0.01; 
Xmax = 5; 
SOL 2; 

p = TaWe[{A, N[Im[M/.NSolve[-Exp[I2TT/a]M^/"X 
== n?-^'"(M+l/m)(M+l), M\][[SOL]] / {tiX)]} , 
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{A, width, Xmax, width}]; 
ListPlot[Abs[p]] 

The constant SOL is a positive integer that indicates 
which of the possible solutions to pick. A value of a not ex- 
pressed as a fraction of integers causes a dramatic increase 
in running time, which otherwise is less than a minute. 

The Monte Carlo approximation of a free stable ran- 
dom matrix A described in Sec. \3\ its use to build a free 
Wishart-Levy matrix C, and the numerical computation 
of the eigenvalue spectrum of the latter including the sta- 
tistical averaging described in Sec. [5] were carried out with 
Matlab 7.5: 

alpha = 3/2; 7- index of Levy stable distribution 

gam =1; 7- scale parameter of Levy stable distribution 

width = .05; Z bin width of eigenvalue histogram 

N = 200; 7o number of time series 

T = 600; 7o points in each time series; must be >= N. 

R =20; 7o rcindom rotations 

S = 36000; 7o number of sampled eigenvalues 

psi = CT*R*gammaCl+alpha) ) ~ (2/alpha) ; 7. normalisation factor 
rho = [] ; 7o set up array of eigenvalues 
iS = 0; 7o initialise normalisation counter 

while CiS < S) 

7a approximation of a free stable matrix 
L = stabrnd(alpha,0,gam,0,T,T) ; 
for iR = 2:R 

[0,U] = qr(rcindn(T,T)) ; 7= is a random orthogonal matrix 
L = L + 0*stabrndCalpha,0,gam,0,T,T)*0' ; 

end 

7o average over covariance matrices 
for i = 1:N:T-N+1 

Mi = L(i : i+K-1 , : ) ; % choose K out of T rows from L 

Ci = Mi*Mi'/psi; 7= normalisation 

rho = [rho eig(Ci)']; 7» collect the eigenvalues 

iS = iS + N; 

if (iS >= S) 
break; 

end 

end 

end 

[histrho Irho] = hist (rho , : width : 100) ; 7o build the histogram 
histrho = histrho/ Clength(rho) *width) 7. normalisation 
'/o Irho contains the abscissa and histrho the ordinate 

On a 2.2 GHz AMD Athlon 64 X2 "Toledo" Dual-Core 
with Fedora Core 7 Linux, all the Monte Carlo calculations 
for Fig. [1] together lasted about 6.6 hours, ranging from 
less than 2 minutes each for a = 1, to = 1 to about 47 
minutes for a ^ 1, to = 1/6. The slow step is the ap- 
proximation of A, i.e. the first for-loop, while the second 
for-loop with the diagonalisation takes from a maximum 
of 2.5% of the total time for a = 1, to = 1 down to 0.25% 
for a ^ 1, to = 1/6. This matches the observation, which 
we made in the range = 10-800 and for the values of 
a, TO, R, S reported in Fig. [l] that the CPU time is ap- 
proximatively proportional to = (N/m)'^ and lower for 
a = 1. In this case, corresponding to the Cauchy distri- 
bution, Eq. ([45]) reduces to X ~ jtancj), which requires 
fewer operations than the general formula. 
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